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ABSTRACT 

We explore how wave-particle interactions affect diffusive shock acceleration (DSA) at astrophysical 
shocks by performing time-dependent kinetic simulations, in which phenomenological models for mag- 
netic field amplification (MFA), Alfvenic drift, thermal leakage injection, Bohm-like diffusion, and a free 
escape boundary are implemented. If the injection fraction of cosmic-ray (CR) particles is £ > 2 x 10 -4 , 
for the shock parameters relevant for young supernova remnants, DSA is efficient enough to develop a 
significant shock precursor due to CR feedback, and magnetic field can be amplified up to a factor of 20 
via CR streaming instability in the upstream region. If scattering centers drift with Alfven speed in the 
amplified magnetic field, the CR energy spectrum can be steepened significantly and the acceleration 
efficiency is reduced. Nonlinear DSA with self-consistent MFA and Alfvenic drift predicts that the 
postshock CR pressure saturates roughly at ~ 10 % of the shock ram pressure for strong shocks with a 
sonic Mach number ranging 20 <; M s <, 100. Since the amplified magnetic field follows the flow modi- 
fication in the precursor, the low energy end of the particle spectrum is softened much more than the 
high energy end. As a result, the concave curvature in the energy spectra does not disappear entirely 
even with the help of Alfvenic drift. For shocks with a moderate Alfven Mach number (Ma < 10), 
the accelerated CR spectrum can become as steep as E~ 2A — E~ 23 , which is more consistent with the 
observed CR spectrum and gamma-ray photon spectrum of several young supernova remnants. 
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1. INTRODUCTION 

Diffusive shock acceleration (DSA) theory explains 
how nonthermal particles are produced through their 
interactions with MHD waves in the converging flows 
across collisionless shocks in astrophysical plasmas 
(Bell 1978; Drury 1983; Blandford & Eichler 1987). 
Theoretical studies have shown that some suprathcr- 
mal particles with velocities large enough to swim 
against the downstream flow can return across the 
shock and stream upstream, and that streaming mo- 
tions of high energy particles against the background 
fluid generate both resonant and nonresonant waves 
upstream of the shock (Bell 1978; Lucek & Bell 2000; 
Bell 2004; Riquelme & Spitkovsky 2009; Rogachevskii 
et al. 2012). Those waves in turn scatter CR particles 
and amplify turbulent magnetic fields in the preshock 
region. These plasma physical processes, i.e., injection 
of suprathermal particles into the CR population, self- 
excitation of MHD waves, and amplification of mag- 
netic fields are all integral parts of DSA (e.g., Malkov 
& Drury 2001). 

Multi-band observations of nonthermal radio to 7- 
ray emissions from supernova remnant (SNR) shocks 
have confirmed the acceleration of CR electrons and 
protons up to - 100 TeV (e.g., Abdo ct al. 2010, 2011; 
Acero et al. 2010; Acciari et al. 2011; Giordano et al. 



2012). Moreover, thin rims of several young SNRs in 
high-resolution X-ray observations indicate the pres- 
ence of downstream magnetic fields as strong as a few 
100/iG, implying efficient magnetic field amplification 
(MFA) at these shocks (e.g., Parizot et al. 2006; Erik- 
sen et al. 2011; Reynolds et al. 2012). 

The most attractive feature of the DSA theory is the 
simple prediction of power-law energy spectra of CRs, 
N(E) cx _E _ ( <T + 2 )/( <T_:L ) (where a is the shock compres- 
sion ratio) in the test particle limit. For strong, adi- 
abatic gas shocks with a = 4, this gives a power-law 
index of 2, which is reasonably close to the observed 
'universal' index of the CR spectra in many environ- 
ments. However, nonlinear treatments of DSA pre- 
dict that at strong shocks there are highly nonlinear 
back-reactions from CRs to the underlying flow, cre- 
ating a shock precursor (e.g., Berezhko & Volk 1997; 
Kang & Jones 2007). So the particles just above the 
injection momentum (p; n j) sample mostly the compres- 
sion across the subshock (a s ), while those near the 
highest momentum (p m ax) experience the greater, to- 
tal compression across the entire shock structure (cr*). 
This leads to the CR energy spectrum that behaves 
as N(E) cx £-(^+2)/(t 3 -i) for p _ p . nj) but flattens 

gradually to N(E) cx eH^ 2 )/^*- 1 ) toward p ~ p max 
(Kang et al. 2009). For example, the power-law index 
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becomes 1.5 for cr t = 7. 

In contrast to such expectations, however, the GeV- 
TeV 7-ray spectra of several young SNRs seem to re- 
quire the proton spectrum as steep as N(E) oc E~ 2 - 3 , 
if the observed 7-ray photons indeed originate from ir° 
decay (Abdo et al. 2010; Giordano et al. 2012). This is 
even softer than the test-particle power-law for strong 
shocks. Moreover, Ave et al. (2009) showed that the 
spectrum of CR nuclei observed at Earth can be fitted 
by a single power law of J(E) oc E~ 2 - 67 below 10 14 
cV. Assuming an energy-dependent propagation path 
length (A oc E~ - 6 ), they suggested that a soft source 
spectrum, N(E) oc E~ a with a ~ 2.3 — 2.4 is preferred 
by the observed data. These observational data ap- 
pear to be inconsistent with flat CR spectra predicted 
by nonlinear DSA model for the SNR origin of Galactic 
CRs. 

It has been suggested that non-linear wave damp- 
ing and wave dissipation due to ion-neutral collisions 
may weaken the stochastic scattering on relevant scales, 
leading to slower acceleration than predicted based on 
the the so-called Bohm-like diffusion, and escape of the 
highest energy particles from the shock (e.g. Ptuskin & 
Zirakashvili 2005; Caprioli et al. 2009). These processes 
may lead to the particle energy spectrum at the highest 
energy end that is much steeper than predicted by non- 
linear DSA. Escape of high energy protons from SNRs 
is an important yet very complex problem that needs 
further investigation (Malkov et al. 2011; Drury 2011). 

Recently some serious efforts have been underway 
to understand at least some of the complex plasma 
processes through Particle-in-Cell (PIC) and hybrid 
plasma simulations (e.g. Riquelme & Spitkovsky 2009; 
Guo et al. 2012; Garate & Spitkovsky 2012). How- 
ever, these types of plasma simulations are too much 
demanding and too expensive to study the full extent of 
the DSA problem. So we do not yet understand them in 
enough detail to make precise quantitative predictions 
for the injection and acceleration rate and efficiency. 
Instead, most of kinetic approaches commonly adopt 
phenomenological models that can emulate more or less 
self-consistently some of those plasma interactions, for 
example, the thermal leakage injection, magnetic field 
amplification, wave-damping and etc (e.g., Berezhko et 
al. 2009; Kang 2010; Ptuskin et al. 2010; Lee et al. 2012; 
Caprioli 2012). 

In our previous studies, we considered DSA of 
CR protons, assuming that magnetic field strength is 
uniform in space and constant in time without self- 
consistent MFA (e.g. Kang & Jones 2007; Kang et al. 
2009). In the present paper, we explore how the fol- 
lowing processes affect the energy spectra of CR pro- 
tons and electrons accelerated at plane astrophysical 
shocks: 1) magnetic field amplified by CR stream- 
ing instability in the precursor 2) drift of scattering 
centers with Alfven speed in the amplified magnetic 
field, and 3) escape of highest energy particles from 
the shock. Toward this end we have performed time- 



dependent numerical simulations, in which DSA of CR 
protons and electrons at strong planar shocks is fol- 
lowed along with electronic synchrotron and inverse 
Compton (IC) losses. Magnetic field amplification due 
to resonant waves generated by CR streaming instabil- 
ity is included through an approximate, analytic model 
suggested by Caprioli (2012). Escape of highest en- 
ergy particles near maximum momentum, p max , is in- 
cluded by implementing a free escape boundary (FEB) 
at a upstream location. As in our previous works (e.g. 
Kang 2010, 2011), a thermal leakage injection model, a 
Bohm-likc diffusion coefficient (n(p) oc p), and a model 
for wave dissipation and heating of the gas are adopted 
as well. 

In the next section we describe the numerical method 
and phenomenological models for the plasma interac- 
tions in DSA theory, and the model parameters for pla- 
nar shocks. The simulation results will be discussed in 
Section 3, followed by a brief summary in Section 4. 

2. DSA MODEL 

2.1 CRASH Code for DSA 

Here we consider the CR acceleration at quasi- 
parallel shocks where the mean background magnetic 
field lines are parallel to the shock normal. So we 
solve the standard gasdynamic equations with CR pro- 
ton pressure terms added in the conservative, Eulerian 
formulation for one dimensional plane-parallel geome- 
try. The basic gasdynamic equations and details of the 
CRASH (Cosmic-Ray Amr SHock) code can be found 
in Kang et al. (2002) and Kang (2011). 

We solve the following diffusion-convection equa- 
tions for the pitch-angle-averaged phase space distri- 
bution function for CR protons, g p = f p p 4 , and for CR 
electron, g e — f e p 4 (Skilling 1975): 
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where y = \n(p/m p c). Here the particle momentum is 
expressed in units of m p c and so the spatial diffusion 
coefficient, n(x,p), has the same form for both protons 
and electrons. The velocity u w represents the effective 
relative motion of scattering centers with respect to the 
bulk flow velocity, u, which will be described in detail 
in section 2.5. 

The cooling term b(p) = —dp/dt takes account for 
electron synchrotron/IC losses, while it is set to be 
b(p) = for protons. Here the synchrotron/IC cool- 
ing constant for electrons is defined as 



b(p) 



4e 4 



9m 4 c 6 



bIp 2 



(2) 
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in cgs units, where e and m e are electron charge and 
mass, respectively. Here B 2 = B 2 + B 2 as the effective 
magnetic field strength for radiative losses including 
the energy density of the ambient radiation field. We 
set B r = 6.5 p,G, including the cosmic background and 
mean Galactic radiation fields (Edmon et al. 2011). 

The dynamical effects of the CR proton pressure are 
included in the DSA simulations, while the CR elec- 
trons are treated as test-particles. In order to include 
the dynamical effects of amplified magnetic field, the 
magnetic pressure, Pb = B 2 /8tt, is added to the mo- 
mentum conservation equation as follows: 



d{pu) d(pu 2 



P g +Pc + P B ) 



dt 



dx 



= 0. 



(3) 



However, our model magnetic field amplification typi- 
cally results in Pb/ pou 2 < 0.01 in the precursor, where 
poll 2 is the shock ram pressure (see Section 2.4). 

2.2 Thermal Leakage Injection 

The injection rate with which suprathermal parti- 
cles are injected into CRs at the subshock depends in 
general upon the shock Mach number, field obliquity 
angle, and strength of Alfvenic turbulence responsi- 
ble for scattering. In thermal leakage injection models 
suprathermal particles well into the exponential tail of 
the postshock Maxwellian distribution leak upstream 
across a quasi-parallel shock (Malkov & Drury 2001; 
Kang et al. 2002). We adopt a simple injection scheme 
in which the particles above an effective injection mo- 
mentum pi„j cross the shock and get injected to the CR 
population: 



Pinj « 1. 17 m p u 2 1 



1.07 



(4) 



where the injection parameter, es = Bq/Bj_, is the 
ratio of the large-scale magnetic field along the shock 
normal, B , to the amplitude of the postshock MHD 
wave turbulence, B± (Kang et al. 2002). With a larger 
value of 6b (i.e., weaker turbulence), p- m j is smaller, 
which results in a higher injection rate. We consider 
e B = 0.23 here. 

We define the injection efficiency as the fraction of 
particles that have entered the shock from far upstream 
and then injected into the CR distribution: 



K e/p ~ 10~ 4 - 10~ 2 (Reynolds 2008; Morlino & Capri- 
oli 2012). Since this ratio is not yet constrained accu- 
rately by plasma physics and we do not consider non- 
thermal emissions from CR particles in this paper, both 
protons and electrons arc injected in the same man- 
ner in our simulations (i.e. basically K e / p = 1). But 
K e / P = 0.1 will be used just for clarity of some figures 
below. 

2.3 Bohm-like Diffusion Model 

It is assumed that CR particles are resonantly scat- 
tered by Alfven waves, which are excited by CR stream- 
ing instability in the upstream region and then ad- 
vected and compressed in the down stream region (Bell 
1978; Lucek & Bell 2000). So in DSA modeling the 
Bohm diffusion model, Kb = (l/3)r g v, is commonly 
used to represent a saturated wave spectrum. We adopt 
a Bohm-like diffusion coefficient that includes a flat- 
tened non-relativistic momentum dependence, 



k(x,p) = K n 



B fl 



P 



B(x) m v c 



(6) 



where n n = m p c 3 /(3eB ) = (3.13 x 10 22 cm 2 s _1 ).Bo S 
and Bq is the magnetic field strength far upstream ex- 
pressed in units of microgauss. The local magnetic field 
strength, B{x), will be described in the next section. 
Hereafter we use the subscripts '0', '1', and '2' to denote 
conditions far upstream of the shock, immediately up- 
stream and downstream of the subshock, respectively. 

2.4 Magnetic Field Amplification 

Since the resonant interactions amplify mainly the 
turbulent magnetic field perpendicular to the shock 
normal in the quasi-linear limit, it was commonly as- 
sumed that the parallel component is Bu w Bq, the 
unperturbed upstream field (Caprioli et al. 2009). In 
a strong MFA case, however, the wave-particle inter- 
action and the CR transport arc not yet understood 
fully. For example, plasma simulations by Riquelme & 
Spitkovsky (2009) showed that both B\\/B and B ± /B Q 
can increase to ~ 10 — 30 via Bell's CR current driven 
instability. Here we follow the prescription for MFA 
that was formulated by Caprioli (2012) based on the as- 
sumption of isotropization of amplified magnetic field. 

In the upstream region (x > x s ), 



J dx J* p °° 4tt/ p (p, x, t)p 2 dp 
n u s t 



(5) 



where n is the particle number density far upstream 
and u s is the shock speed. 

Since postshock thermal electrons need to be pre- 
accelerated before they can be injected into Fermi pro- 
cess, it is expected that electrons are injected at the 
shock with a much smaller injection rate, i.e., the 
CR electron-to-proton ratio is estimated to be small, 



B{xf 



1 + (1 - uh) 



25 



-M 



A,0' 



(1 - U(xfl *) 
U(x) 3 / 2 



U2 



(7) 



where Ma,o = u s /va.o is the Alfven Mach number for 
the far upstream Alfven speed, va,o = Bo/y/Anpo, and 
U(x) = [u s — u(x)]/u s is the flow speed in the shock 
rest frame normalized by the shock speed. The factor 
(1 — u>h) is introduced to take account of the loss of 
magnetic energy due to wave dissipation, which will be 
discussed in Section 2.5. Obviously, ljh = means no 
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dissipation, while ujh = 1 means complete dissipation 
of waves (i.e., no MFA). Here loh — 0.1 will be consid- 
ered as a fiducial case, since we are interested in the 
case where the effects of MFA and ensuing wave drift 
are the greatest. 

This MFA model predicts no amplification in the 
test-particle regime, where the flow structure is not 
modified (i.e., U(x) = 1). In the case of "moder- 
ately modified" shocks, for example, if U\ k> 0.8 and 
ojh =0.1, the amplified magnetic field strength scales 
as B 1 /B Q w 0.11M A)0 . So for M Afi w 150, the 
preshock amplification factor could become Bi/B w 
17. On the other hand, the ratio of the magnetic pres- 
sure to the shock ram pressure becomes Pb,i/po u I = 

(2/25)(l-?7 1 V4 ) 2 /?7 1 3/2 w 6.6xl0~ 3 . So we expect that 
even the amplified field is not dynamically important 
in the precursor. 

The magnetic filed strength immediately upstream 
of the subshock, B\, is estimated by Equation (7) and 
assumed to be completely turbulent. Moreover, assum- 
ing that the two perpendicular components are simply 
compressed across the subshock, the immediate post- 
shock field strength can be estimated by 



B 2 /B 1 = v / l/3 + 2/3( (02 / /9l ) 2 . 



(8) 



So for the case with p 2 /pi ~ 4.2, B 2 jB\ k> 3.5. Then 
we assume in the downstream region the field strength 
scales with the gas density: 



B(x) = B 2 ■ [p(x)/p 2 ] 



(9) 



We note the MFA model described in Equations (7)- 
(9) is used also for the diffusion coefficient model given 
in Equation (6). Hence the maximum momentum p max 
is controlled by the degree of MFA as well. 

2.5 Alfvenic Drift 

The resonant waves generated by CR streaming in- 
stability will drift with respect to the underlying flow 
and also transfer energy to the gas through dissipation 
(e.g. Skilling 1975; Jones 1993). These two effects in- 
fluence the accelerated particle spectrum and the DSA 
efficiency as follows. The scattering by Alfven waves 
tends to isotropize the CR distribution in the wave 
frame rather than in the gas frame (Bell 1978), which 
reduces the velocity difference between upstream and 
downstream scattering centers, compared to that of the 
bulk flow. The resulting CR spectrum becomes softer 
than estimated without considering the wave drift. 

The mean drift speed of scattering centers is com- 
monly set to be the Alfven speed, i.e., u wA (x) = +va = 
-\-Bu/ y/A-Kp, pointing away from the shock, where Bu is 
the local magnetic field strength parallel to the shock 
normal. As described in Equation (7) here we assume 
that both B\\ and B± are amplified and isotropized, so 
scattering centers drift with Alfven speed in the local 
amplified magnetic field. In order to take account of 



the uncertainty regarding this issue, we model the local 
Aflven speed as 



v A (x) 



Bg + (B(x) - B )f A 

\/^Kp{x) 



(10) 



where the parameter f A is a free parameter (Zi- 
rakashvili & Ptuskin 2008; Lee et al. 2012). If scat- 
tering centers drift along the amplified field (f A = 1), 
the Alfvenic drift will have the maximum effects. Here 
we will consider the models with f A = 0.5 — 1.0 (see 
Table 1). 

In the postshock region the Alfvenic turbulence is 
probably relatively balanced, so the wave drift can be 
ignored, that is, u w , 2 ~ (Jones 1993). On the other 
hand, if the scattering centers drift away from the shock 
in both upstream and downstream regions, the accel- 
erated particle spectrum could be softened drastically 
(e.g. Zirakashvili & Ptuskin 2008). We will consider 
one model (H2d) in which u w , 2 ss — va is adopted in 
the downstream of the shock (see Table 1). 

As mentioned in the Introduction, the CR spectrum 
develops a concave curvature when the preshock flow is 
modified by the CR pressure. If we include the Alfvenic 
drift only in the upstream flow, the slope of the mo- 
mentum distribution function, q = —dlnf/dlnp, can 
be estimated as 



3(iti - u w> i) 



2,a s {l-M A \) 



(ui-u w ,i)-u 2 <7 S (1 - M A \) - 1 



for p <~ pi n j , and 



Qt 



3(u - Uwfl) 



3a t (l-M^ ) 



(U -U w fl)-U 2 <7t(l - M A ] ) - 1 



(11) 



(12) 



for p ~ Pmax- Here Ma,i = Ui/v A ,i is Alfvenic Mach 
number immediately upstream of the subshock. As can 
be seen in these equations, a significant steepening will 
occur only if M A £ 10 (Caprioli 2012). 

According to the MFA prescription given in Equa- 
tion (7), the amplification factor depends on the pre- 
cursor modification, that is, the ratio B(x)/B is unity 
far upstream and increases through the precursor to- 
ward the subshock. So the Afvcnic drift speed is high- 
est immediately upstream of the subshock, while it is 
the same as the unperturbed Alfven speed, va,o at the 
far upstream region {Ma,i <C M a ,o)- Thus the Alfvenic 
drift is expected to steepen preferentially the lower en- 
ergy end of the CR spectrum, since the lowest energy 
particles diffuse mostly near the subshock. For the 
highest energy particles, which diffuse over the distance 
of ~ K(p P ,max)/ws, however, the Alfevnic drift does not 
steepenthe CR spectrum significantly, if Ma,o 3> 1- 

2.6 Wave Dissipation and Particle Escape 

As discussed in the Introduction, non-linear wave 
damping and dissipation due to ion-neutral collisions 
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Table 1. Model Paramcters a 



Modcl b 


U s 


n H (ISM) 
(cm -3 ) 


To 


M s 


M A .o 


fA c 


ojh d 


u w ,i 


U W} 2 




km s -1 


(K) 














Wla 


3 x 10 3 


1.0 


4.0 x 10 4 


100 


164. 


1.0 


0.1 


+VA 





Wlb 


3 x 10 3 


1.0 


4.0 x 10 4 


100 


164. 


1.0 


0.1 








Hla 


3 x 10 3 


1.0 


10 6 


20 


164. 


1.0 


0.1 


+VA 





Hlb 


3 x 10 3 


1.0 


10 6 


20 


164. 


1.0 


0.1 








Hlc 


3 x 10 3 


1.0 


10 6 


20 


164. 


0.5 


0.5 


+VA 





H2a 


3 x 10 3 


0.01 


10 6 


20 


16.4 


1.0 


0.1 


+VA 





H2b 


3 x 10 3 


0.01 


10 6 


20 


16.4 


1.0 


0.1 








H2d 


3 x 10 3 


0.01 


10 6 


20 


16.4 


1.0 


0.1 


+VA 


-VA 


H3a 


10 3 


0.01 


10 6 


6.67 


5.46 


1.0 


0.1 


+VA 





H3b 


10 3 


0.01 


10 6 


6.67 


5.46 


1.0 


0.1 








H4a 


4.5 x 10 3 


0.01 


10 6 


30 


24.6 


1.0 


0.1 


+VA 





H4b 


4.5 x 10 3 


0.01 


10 6 


30 


24.6 


1.0 


0.1 









a For all the models the background magnetic field is Bo — 5 /iG and the injection parameter is e# = 0.23. 

b 'W and 'H' stands for the warm and hot phase of the ISM, respectively. 

c See (9) for the Alfven parameter. 

d See Equation (6) for the wave dissipation parameter. 



may weaken the stochastic scattering, leading to slower 
acceleration and escape of highest energy particles from 
the shock. These processes are not understood quan- 
titatively well, so we adopt a simple model in which 
waves are dissipated locally as heat in the precursor. 
Then gas heating term in the upstream region is pre- 
scribed as 



dP c 

W(x,t) = -lu h ■ v A (x)-^-, 



(13) 



where P c is the CR pressure (Jones 1993). The pa- 
rameter ujh is introduced to control the degree of wave 
dissipation and a fiducial value of ujh = 0.1 is assumed. 
As shown previously in SNR simulations (e.g. Bcrezhko 
& Volk 1997; Kang & Jones 2006), this precursor heat- 
ing reduces the subshock Mach number thereby reduc- 
ing the DSA efficiency. For larger values of ujh, the 
magnetic field amplification is suppressed (see Equa- 
tion (7)), which also reduces the maximum momentum 
of protons and so the DSA efficiency. 

In addition, we implement a free escape boundary 
(FEB) at a upstream location by setting J(xfeb,p) = 
at xfeb = O.lRg = 0.3pc (here the shock is located 
at x s = 0). This FEB condition can mimic the es- 
cape of the highest energy particles with the diffusion 
length, k(p)/u s J> xfeb- For typical supernova rem- 
nant shocks, this FEB leads to the size-limited maxi- 
mum momentum, 



4.4xl0 4 (-^)( 



5^G"3000 km s-^O.Spc 



As can be seen in Section 3, the CR proton spectrum 
and the shock structure approach to time-asymptotic 
states, if this FEB is employed (Kang et al. 2009). 



On the other hand, the maximum electron momen- 
tum can be estimated by 



« 2.8 x 10 4 



W 2 ( 

\30nGj V 



3000 km s- 



which is derived from the equilibrium condition that 
the DSA momentum gains per cycle are equal to the 
synchrotron/IC losses per cycle (Kang 2011). The elec- 
tron spectrum at the shock position, f e {x s ,p), cuts off 
exponentially at <~ p ,max- 

On the other hand, the postshock electron spec- 
trum cuts off at a progressively lower momentum down- 
stream from the shock due to the energy losses. That 
results in the steepening of the volume integrated elec- 
tron energy spectrum, F e (p) = j f e (x,p)dx, by one 
power of the momentum (Kang et al. 2012). At the 
shock age t, the break momentum can be estimated 
from the condition t — p/b(p): 



^^«1.3xl0 3/ ' 



10 3 yr 



B. 



e,2 



100 ^G 



, (16) 



which depends only on the postshock magnetic field 
strength and the shock age (Kang 2011). 

2.7 Planar Shock Parameters 

We consider planar shocks with u s — 1000 — 
4500 kms -1 , propagating into a uniform ISM mag- 
netized with Bo = 5 ^G. The model parameters are 
summarized in Table 1. Previous studies have shown 
that the shock sonic Mach number is one of the key pa- 
rameter governing the evolution and the DSA efficiency 
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■*| — i— rn — | — i — i — i — | — rn — i - 1 — i — i — i — | — i — r 

Hlb / / ; 



-Hla 



J i i i I i i i L 



J I I ' I I L 



-0.6 -0.4 -0.2 0.2 
x(pc) 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 


i i i 1 i i 


- /- / / , /! 




-M A0 =164 1 




ji H =lcm~ 3 ; 




-u s =3000km/s 




"T =10 6 K 
i 1 i i i 1 i i i 1 i i i 1 


i i 1 i i 



1 1 II 1 II 1 1 II II 1 1 1 

rt/t n =o.i, o.5, i.o ~ 


■ i 1 i i i r ■ 


i i 1 i 1 1 





-0.6 -0.4 -0.2 0.2 

x(pc) 
log(G p /t), log(G e /t) 



i — i — i — I — i — i — i — I — i — i — i — I — i — r 



0.6 



-0.4 -0.2 
x(pc) 



0.2 




2 4 
log(p/m p c) 



Fig. 1. — Time evolution of the magnetic field strength, CR pressure, gas density, and volume integrated distribution 
functions of protons (G p ) and electrons (G e ) for Hla (solid lines) and Hlb (dotted lines) models at t/t n — 0.1, 2.5 and 5. 
See Table 1 for other model parameters and normalization constants. In the bottom right panel the upper curves are for 
the proton spectra, while the lower curves are for the electron spectra. Note that both G p /t and G e /t are given in arbitrary 
units and K e / p = 0.1 is adopted here for clarity. 



(e.g. Kang & Jones 2007; Kang et al. 2009), so here two 
phases of the ISM are considered: the warm phase with 
T = 4 x 10 4 K (W models), and the hot phase with 
To = 10 6 K (H model). The sonic Mach number of each 
model is given as M s = 20(T /10 6 K) _1 / 2 m 3 ooo, where 
W3000 = Ms/3000 km s -1 . Two values of the gas den- 
sity uh = 0.01 cm -3 and 1 cm -3 are considered. The 
upstream Alfven speed is then va,o — Bq / ^/Airpo = 

(18.3 km s~ 1 )n H 1 ^ 2 , so the Alfvenic Mach number is 
M A ,o = u s /v A ,a = 164^rT7fu 3 ooo- 

We consider Wla, Hla, H2a, H3a and H4a models 
as fiducial cases with canonical values of model param- 
eters: Ja — 1.0 and u>h = 0.1. In models Hlb, H2b, 
H3b and Hlb, Alfvenic drift is turned off (u w> i = 0) 
for comparison. But we note that these models are 
not self-consistent with our MFA model, which assumes 
that Alfven waves propagate along the amplified mag- 
netic field. In Hlc model, MFA is reduced by setting 
Ja — 0.5 and wh = 0.5 . Model H2d is chosen to see 
the effects of Alfvenic drift in the postshock region. 

The physical quantities are normalized in the nu- 



merical code and in the plots below by the follow- 
ing characteristic values: u n = u s , x n = R s =3 pc, 
t n = x n /u n = (978 yr)M 300 o, K n = u n x n , p n = 
(2.34 x 10~ 24 g cm" 3 ) • n H , and P n = p n u 2 n = (2.11 x 
10" 7 erg cm" 3 ) • n H u\ mQ . 

3. DSA SIMULATION RESULTS 

Figures 1-3 show the spatial profiles of the model 
magnetic field, CR pressure, gas density, and the vol- 
ume integrated distribution functions of protons (G p = 
J g p (x,p)dx) and electrons (G e = J g e (x,p)dx) for Hla 
and Hlb (M Afi = 164), H2a and H2b {M A ,o = 16.4), 
and H3a and Hb (Ma,o = 5.46) models, respectively. 
In these simulations, the highest level of refinement is 
Jg.max = 8 and the factor of each refinement is two, 
so the ratio of the finest grid spacing to the base grid 
spacing is Axs/Axo = 1/256 (Kang et al. 2002). Since 
Figures 1-3 show the flow structure on the base grid, 
the precursor profile may appear to be resolved rather 
poorly here. More accurate values of the precursor den- 
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Same as Figure 1 except that H2a (solid lines) and H2b (dotted lines) are shown. 



sity compression (p±) and magnetic field amplification 
(Bi ) can be found in Figure 4 below. Also note that the 
FEB is located at xfeb = 0.3pc and we use K e / p = 0.1 
here in order to show the proton and electron spectra 
together. 

These figures demonstrate that 1) the shock struc- 
ture reaches the time-asymptotic state and evolves in 
a self-similar fashion for t/t n J> 0.1 (Kang & Jones 
2007). 2) Also the proton spectrum approaches the 
steady state for t/t n > 0.5 due to the FEB, while 
the electron spectrum continues to cool down in the 
downstream region. 3) Magnetic field is amplified by a 
greater factor for a higher Ma,o- 4) Alfvenic drift steep- 
ens the CR spectrum, and reduces the CR acceleration 
efficiency and flow modification by CR feedback, result- 
ing in lesser MFA. 5) At low energies the CR spectra 
are much steeper than the test-particle power-law due 
to the velocity profile and magnetic field structure in 
the precursor. 

In Hla model with Alfvenic drift (solid lines), the 
gas density immediately upstream of the subshock in- 
creases to pi/po ~ 1.2, while the total compression ra- 
tio becomes P2/P0 ~ 4.5. So the flow structure is mod- 
erately modified by the CR pressure feedback: U\ ~ 0.9 
and P c ,2/ ' Pqu 2 s ~ 0.12. Since the Alfven Mach number 



is high (Ma,o = 164), the self-amplified magnetic field 
strength, based on the model in equation (7), increases 
to B\ w 100/jG, which results in the immediate post- 
shock fields of B2 ~ 350/xG. Compared to the model 
without Alfvenic drift (Hlb, dotted lines), the CR dis- 
tribution functions are softer. Although Alfvenic drift 
steepens the distribution function in Hla model, G p (p) 
still exhibits a significant concave curvature and it is 
slightly flatter than the test-particle power-law (E~ 2 ) 
at the highest energy end. This is because Ma,o = 164 
is too large to induce the significant enough softening 
for 

P Pp,max (see Equation (12)). Note that Pp^max is 
lower in Hla model than that in Hlb model, because 
of weaker MFA. 

The structures of the integrated electron spectra are 
complex for p <; p e ,max- As shown in Equation (16), the 
volume integrated electron energy spectrum steepens 
by one power of the momentum due to radiative cool- 
ing. One can see that the break momentum, phr(t), 
shifts to lower momenta in time. The peak near p e ,max 
comes from the electron population in the upstream 
region, which cools much less efficiently due to weaker 
magnetic field there (Edmon et al. 2011). Since MFA is 
much stronger in Hlb model, compared to Hla model, 
the electron spectrum cools down to lower momentum 
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in the downstream region. 

Comparing H2a/b in Figure 2 with Hla/b in Fig- 
ure 2, one can see that the degree of shock modifica- 
tion is similar in the these models. Because of a lower 
gas density (njj — 0.01 cm~ 3 ) in H2a/b models, the 
Alfven Mach number is smaller (Ma,o — 16.4) and so 
MFA is much less efficient, compared to Hla/b models. 
In H2a model the amplified preshock field increases to 
only B\ « 10 pG, while the postshock field reaches 
.E?2 ~ 35 /iG. Because of much weaker magnetic field, 
compared to Hla/b model, the electron spectra are af- 
fected much less by radiative cooling. 

Since H3a/b models in Figure 3 have a lower sonic 
Mach Number (M s = 6.7), the flow structures are al- 
most test-particle like with B\ rs Bq, pil Pq ~ 4, and 
Pel I Po u s ~ 0.05 — 0.13. So the CR acceleration, flow 
modification, and MFA are all less efficient, compared 
to Hla/b and H2a/b models. In H3a model, especially, 
the CR spectra are as steep as E~ 2A — E~ 23 and elec- 
trons do not suffer significant cooling. 

Figure 4 shows how various shock properties change 
in time for different models: the CR injection frac- 
tion, postshock CR pressure, density compression fac- 
tors and magnetic field strengths. As discussed above, 
the magnetic field amplification is more efficient in the 



models with higher M A ,o- B 1 /B ss 20 for M A ,o = 164 
(Hla, Hlc, Wla models), B 1 /B » 3 for M A = 24.6 
(H4a), Bx/Bq ps 2 for M A0 = 16.4 (H2a), Bi/B ~ 1 
for M A ,o = 5.46 (H3a). 

According to previous studies, nonlinear DSA with- 
out self-consistent MFA and Alfvenic drift predicts 
that the DSA efficiency depends strongly on the sonic 
Mach number M s and the CR pressure asymptotes to 
P c ,2/Pqu 2 ~ 0.5 for M s > 20. (Kang & Jones 2007; 
Kang et al. 2009). However, Figure 4 shows that in the 
models with MFA and Alfvenic drift the CR accelera- 
tion and MFA are reduced in such a manner that the 
DSA efficiency saturates roughly at P c ,2/po u l ~ 0.1 for 
20 <; M s <; 100. We can see that models with a wide 
range of sonic Mach number, i.e. Wla(M s = 100), H4a 
(M s = 30), H2a and Hla (M s = 20), all have similar 
results: P2/P0 ~ 4.5, and P c .il 'pou 2 ~ 0.1. 

Figures 5 and 6 show the volume-integrated dis- 
tribution function, G p (p)/(noU s t), for protons and, 
G e (p)/(nou s t), for electrons, respectively, for different 
models. Again the proton spectrum approaches to the 
steady state for t/t n ;> 0.5, when p Pimax (t) satisfies 
the condition, K(p P;max )/ws ~ xfeb- We note that 
for Wlb model only the curve at t/t n = 0.1 is shown, 
because the simulation was terminated when the sub- 
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Fig. 4. — Time evolution of the injection efficiency, £, postshock CR pressure, P c ,2, the gas density pi (P2) immediately 
upstream (downstream) of the subshock, and magnetic field strengths, Bi and B2, in different models: Hla (black solid 
lines), Hlb (red dotted), Hlc (blue dashed), Wla (green dot-dashed) in the left column, and H2a (black solid lines), H2b 
(red dotted), H3a (blue dashed), H4a (green dot-dashed) in the right column. Note Hlb and H2b models are shown for 
comparison, but B\ and B2 for those models are not included in the bottom panels. 



shock disappears because of very efficient DSA. These 
figures demonstrate that the CR spectra is steepened 
by Alfvenic drift, especially at lower energies, and that 
the degree of softening is greater for smaller Ma,o- In 
H2d model, in which the downstream drift is included 
(uw,2 = —va) in addition to the upstream drift, the 
CR spectra are steepened drastically. 

In the volume-integrated electron spectrum, the low- 
energy break corresponds to the momentum at which 
the electronic synchrotron/IC loss time equals the 
shock age. In the models with stronger magnetic field 
(e.g., Hla and Wla models), this spectral break occurs 
at a lower p c ,br, and the separate peak around p e ,max 
composed of the upstream population becomes more 
prominent. 



4. SUMMARY 

Using the kinetic simulations of diffusive shock accel- 
eration at planar shocks, we have calculated the time- 
dependent evolution of the CR proton and electron 
spectra for the shock parameters relevant for typical 
young supernova remnants. In order to explore how 
various wave-particle interactions affect the DSA pro- 
cess, we adopted the following phenomenological mod- 
els: 1) magnetic field amplification (MFA) induced by 
CR streaming instability in the precursor, 2) drift of 
scattering centers with Alfven speed in the amplified 
magnetic field, 3) particle injection at the subshock via 
thermal leakage injection, 4) Bohm-like diffusion coeffi- 
cient, 5) wave dissipation and heating of the gas, and 6) 
escape of highest energy particles through a free escape 
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Fig. 5. — Volume integrated distribution function of CR protons for different models: Hla, Wla, H2a, Hlc, H3a, H4a, 
(black solid lines), Hlb, Wlb, H2b, H2d, H3b, H4b (red dotted lines). For Wlb model, only the curve for t/t„ = 0.1 is 
shown, because the simulation was terminated afterward. The curves for H3a/b are multiplied by a factor of 10 to show 
them in the same scale as other models. See Table 1 for model parameters. 



boundary. 

The MFA model assumes that the amplified mag- 
netic field is isotropized by a variety of turbulent pro- 
cesses and so the Alfven speed is determined by the lo- 
cal amplified magnetic field rather than the background 
field (Caprioli 2012). This model predicts the mag- 
netic field amplification factor scales with the upstream 
Alfvenic Mach number as B\jB^ oc Ma,o, and also in- 
creases with the strength of the shock precursor (see 
Equation (7)). 

Moreover, we assume that self-generated MHD waves 
drift away from the shock with respect to the back- 
ground flow, leading to smaller velocity jumps that 
particles experience scattering across the shock. The 
ensuing CR distribution function becomes steeper than 
that calculated without Alfvenic drift, so the CR injec- 
tion/acceleration efficiencies and the flow modification 
due to CR feed back are reduced. 

The expected power-law slope depends on the Alfvenic 
Mach number as given in Equations (11)-(12). With 
our MFA model that depends on the precursor modifi- 
cation, the upstream Alfvenic drift affects lower energy 



particles more strongly, steepening the low energy end 
of the spectrum more than the high energy end. Hence, 
for Ma,o > 10, the CR spectra still retain the concave 

curvature and they can be slightly flatter than E~ 2 at 
the high energy end. For weaker shocks with M s = 6.7 
and Ma,o = 5-5 (H3a model), on the other hand, the 
Alfvenic drift effects are more substantial, so the energy 
spectrum becomes as steep as N(E) oc E~ 21 — E~ 23 . 

We can explain how MFA and Alfvenic drift regu- 
late the DSA as follows. As CR particles stream up- 
stream of the shock, magnetic field is amplified and 
Alfven speed in the local B(x) increases in the precur- 
sor. Then scattering centers drift with enhanced va, 
the CR spectrum is steepened and the CR acceleration 
efficiency is reduced, which in turn restrict the growth 
of the precursor (see also Caprioli 2012). So the flow 
modification due to the CR pressure is only moderate 
with P2/P0 ~ 4.5. As a result, the DSA efficiency sat- 
urates roughly at P c ^/p Q u 2 ~ 0.1 for 20 <; M s < 100. 
For M s = 20 shocks with u s = 3000 km s _1 , for exam- 
ple, in the models with Alfvenic drift (Hla and H2a), 
the CR injection fraction is reduced from ^~2x 10~ 3 
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Fig. 6. — Same as Figure 5 except that the volume integrated distribution function of CR electrons are shown. 



to - 2 x 1CT 4 , while the CR pressure decreases from 
Pc,2/pqu 2 s ~ 0.25 to ~ 0.12, compared to the model 
without Alfvenic drift (Hlb and H2b) (see Figure 4). 

This study demonstrates that detailed nonlinear 
treatments of wave-particle interactions govern the CR 
injection/acceleration efficiencies and the spectra of CR 
protons and electrons. Thus it is crucial to understand 
in a quantitative way how plasma interactions amplify 
magnetic field and affect the transportation of waves 
in the shock precursor through detailed plasma sim- 
ulations such as PIC and hybrid simulations. More- 
over, the time-dependent behaviors of self-amplified 
magnetic field and CR injection as well as particle es- 
cape will determine the spectra of the highest energy 
particles accelerated at astrophysical shocks. We will 
present elsewhere the results from more comprehensive 
DSA simulations for a wide range of sonic and Alfven 
Mach numbers. 
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